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Disclosed are procedures for converting x-ray 
cone beam data (line integrals through an 
object) to Radon data (planar integrals) for 3D 
CT Image reconstruction by inverse Radon 
transformation. The radial derivative of each 
planar integral is detemnined by integrating to 
detenmine weighted line integrals along each of 
a pair of lines on a normalized detector plane, 
which lines are defined as intersections with the 
normalized detector plane of a corresponding 
pair of integration planes sharing a rotation axis 
and rotated with respect to each other by a 
rotation angle 8p, and then dividing the differ- 
ence between the weighted line integrals by the 
rotation angle 8^. The procedure can be em- 
ployed to convert the cone beam date to values 
representing planar integrals on any arbitrary 
set of planes in Radon space, such as a set of 
coaxial vertical planes. 
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Cross-Reference to Related Applications 

The invention disclosed and clainned herein is related to the subject matter of the co-filed patent appli- 
cations, the entire disclosures of which are hereby expressly incorporated herein by reference; based on: 
5 U.S. Serial No.07/631 ,818, filed 21 December 1 990. by Kwok C. Tarn, entitled "PARALLEL PROCES- 

SING METHOD AND APPARATUS FOR RECONSTRUCTING A THREE-DIMENSIONAL COMPUTERIZED 
TOMOGRAPHY (CT) IMAGE OF AN OBJECT FROM CONE BEAM PROJECTION DATA OR FROM PLANAR 
INTEGRALS" [RD-19564]; and 

U.S. Serial No.07/631,514 . filed 21 December 1990, by Kwok C. Tarn, entitled "METHOD AND 
10 APPARATUS FOR RECONSTRUCTING A THREE-DIMENSIONAL COMPUTERIZED TOMOGRAPHY (CT) 
IMAGE OF AN OBJECT FROM INCOMPLETE CONE BEAM PROJECTION DATA" [RD-19695]. 

Background of the Invention 

15 The present invention relates generally to three-dimensional (3D) computerized tomography (CT) and, 
more particularly, to methods and apparatus for converting x-ray cone beam data to planar integrals for 3D 
image reconstruction through inverse Radon transfonmation. 

In conventional computerized tomography for both medical and industrial applications, an x-ray fan beam 
and a linear an-ay detector are employed. Two-dimensional (2D) Imaging is achieved. While the data set is com- 

20 plete and image quality is conrespondingly high, only a single slice of an object is imaged at a time. When a 
3D image is required, a "stack of slices" approach is employed. Acquiring a 3D data set a 2D slice at a time Is 
inherently tedious and time-consuming. Moreover, in medical applications, motion artifacts occur because adja- 
cent slices are not imaged simultaneously. Also, dose utilization is less than optimal, because the distance be- 
tween slices is typically less than the x-ray collimator aperture, resulting in double exposure to many parts of 

25 the body. 

A more recent approach, based on what Is called cone beam geometry, employs a two-dimensional array 
detector Instead of a linear array detector, and a cone beam x-ray source instead of a fan beam x-ray source. 
At any instant the entire object is inradiated by a cone beam x-ray source, and therefore cone beam scanning 
is much faster than sllce-by-slice scanning using a fan beam or a parallel beam. Also, since each "point" in the 

30 object is viewed by the x-rays In 3D rather than in 2D, much higher contrast can be achieved than Is possible 
with conventional 2D x-ray CT. To acquire cone beam projection data, an object Is scanned, preferably over a 
360*' angular range, either by moving the x-ray source in an appropriate scanning trajectory, for example, a 
circular trajectory around the object, while keeping the 2D anray detector fixed with reference to the source, or 
by rotating the object while the source and detector remain stationary. In either case, it is relative movement 

35 between the source and object which effects scanning. 

Most image reconstruction procedures in x-ray CT are based on the Radon inversion process, in which 
the image of an object is reconstructed from the totality of the Radon transform of the object The Radon trans- 
form of a 2D object consists of integrals of the object density on lines intersecting the object The Radon trans- 
fomn of a 3D object consists of planar integrals. The cone beam data, however, are not directly compatible with 

40 image reconstmction through inverse Radon transformation, which requires the use of planar integrals of the 
object as input Consequently, Image reconstruction by inversion from cone beam scanning data generally com- 
prises two steps: (1 ) convert the cone beam data to planar integrals, and (2) perfonm an inverse Radon transfomi 
on the planar integrals to obtain the image. The present invention is primarily directed to efficient nnethods and 
apparatus for converting x-ray cone beam data to planar integrals, or values representing planar integrals, on 

45 a set of arbitrary planes in Radon space. The above-Incorporated application Serial No. 

[RD-19564] discloses a two-step method for performing an inverse Radon transfonm starting with planar integ- 
rals on a set of coaxial vertical planes in Radon space. Thus the invention disclosed herein may be employed 
to convert x-ray cone beam data to values representing planar integrals on a set of coaxial vertical planes in 
Radon space, and the invention of application Serial No. [RD-19564] may be employed 

50 to perform the inverse Radon transformation portion of the 3D image reconstruction. 

One method for converting cone beam data to planar integrals is disclosed In Gerald N. Minerbo, "Con- 
volutional Reconstruction from Cone-Beam Projection Data", IEEE Trans. Nucl. Sci., Vol. NS-26, No. 2, pp. 
2682-2684 (April 1979). Unfortunately, as is discussed, for example, in LA Feldkamp, L.C. Davis, and J.W. 
Kress, "Practical Cone-Beam Algorithm", J. Opt Soc. Am. A., Vol. 1 , No. 6, pp. 612-819 (June 1 984), the deri- 

55 vation In Minerbo contains an error which cannot easily be rectified and which renders the result invalid. 

In Bruce D. Smith, "image Reconstmction from Cone-Beam Projections: Necessary and Sufficient Condi- 
tions and Reconstruction Methods", IEEE Trans. Med. Imag., Vol. MI-44, pp. 1425 (March 1985), there is dis- 
closed a method for converting from cone beam data the one-dimensional convolution of the planar Integrals 
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with the Horn's kernel. Since the convolution mixes together the planar integrals on all the planes, the compu- — ■ - 
tation of one point of the convolved result requires all the data on the detector at one view angle. Thus the task 
is very computationally intensive. 

In P. Grangeat "Analysis of A 3D Imaging System by Reconstruction from X Radiographics in Conical 

5 Geometry" ("Analyse d'un System D-lmagerie 3D par Reconstruction a partir de Radiographies X en Geometrie 
conique"), Ph.D. Thesis, National College of Telecommunications (l-Ecole Nationale Superieure des Telecom- 
. municatlons). France (1 987), a technique is disclosed for computing the derivative of the planar integrals from 
cone beam data. The computed data points, however, reside on a set of great circles on a spherical shell in 
Radon space. These great circles in general do not fall on any arbitrary set of planes In Radon spaces, and do 

10 not fall on a set of coaxial vertical planes in Radon space. Thus they are not suitable for input to inverse Radon 
transfonmation. It would require an extensive effort in three-dimensional interpolation to get the data on the ver- 
tical planes to be used in Inverse Radon transformation, and furthenmore interpolation would introduce errors 
into the data. 

IS Summary of the Invention 

Accordingly, it is an object of the invention to provide methods and apparatus for converting x-ray cone 
beam data to values representing planar integrals on any arbitrary set of planes in Radon space for 3D image 
reconstruction through inverse Radon transformation. 
20 It is a more specific object of the invention to provide methods and apparatus for converting x-ray cone 
beam data to values representing planar integrals on a set of coaxial planes In Radon space for 3D Image recon- 
struction through inverse Radon transfonmation. 

It is another object of the invention to provide such methods and apparatus which require, for each view 
angle, only two closely spaced lines of data on the detector to compute a value representing a planar integral, 
25 in contrast to prior art approaches requiring all data on the detector. 

It is another object of the invention to provide methods and apparatus which are exact and do not require 
interpolation for converting x-ray cone beam data to values representing planar integrals on a set of coaxial 
planes in Radon space, or any arbitrary set of planes in Radon space, for 3D image reconstruction through 
inverse Radon transformation. 
30 It is yet another object of the invention to provide methods and apparatus which minimize the amount of 
computation required to convert x-ray cone beam data to values representing planar integrals on a set of coaxial 
planes in Radon space, or any arbitrary set of planes in Radon space, for 3D image reconstruction through 
Inverse Radon transfomnation. 

In accordance with the invention, there is provided a method for reconstructing a 3D Image of an object 
35 from cone beam projection data, where the cone beam projection data is in the form of line integrals through 
the object organized, for each of a plurality of x-ray source positions S/, as a 20 data set on a nonmalized detec- 
tor plane containing an origin and perpendicular to a line from each particular source position Si to the origin. 
The method Includes the two overall steps of determining values representing planar Integrals on a set of planes 
In Radon space, and then perfonming an inverse Radon transform on the values representing planar integrals 
40 on the set of planes to reconstruct an image of the object. In more particular embodiments, the planes (j>j com- 
prise a set of coaxial planes containing a reference axis intersecting the origin. 

A significant aspect of the invention is the determination of the value of a planar integral or Radon datum 
(actually the radial derivative of the Radon datum) at a particular point in Radon space by Integrating to deter- 
mine weighted line integrals along each of a pair of lines on the nonmalized detector plane, which lines are 
45 defined as intersections with the nonmalized detector plane of a conresponding pair of integration planes sharing 
a rotation axis and rotated with respect to each other by a rotation angle and then dividing the difference 
between the weighted line integrals by the rotation angle 6p. 

Specifically, the step of determining values representing planar integrals on a set of planes <|>j in radon space 
includes the nested steps of, for each of the source positions S/, 
so defining in Radon space a corresponding spherical shell on which Radon data can be determined, 

intersections of the planes (|>j with the spherical shell corresponding to the particular source position S/ defining 
a set of circles D^on the spherical shell, and 
for each of the circles D/y, 

defining a rotation axis as a line through the particular source position S/, intersecting the 
55 particular circle Dy, and orthogonal to the plane of the particular circle Dy, 

defining a set of coaxial Integration planes O/y^ each of the integration planes Q^^ containing 
the particular rotation axis and intersecting the particular circle to define the location of a Radon datum point 
%, and the integration planes intersecting the normalized detector plane on respective lines L^, and 
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— for each of the lines Lfj^ on the normalized detector plane, 

rotating the corresponding integration plane by a small rotation angle 6p to define a 
plane Q|^' Intersecting the normalized detector plane on a corresponding line Lgit\ 

Integrating along the lines and Lgj; to determine respective weighted line integrals J^, 

5 and Jgic\ and 

dividing the difference between the weighted line Integrals and Jgk' by the rotation angle 
. 6p to yield the radial derivative of the Radon datum at the particular point Rgif. 

Similarly, apparatus in accordance with the invention for reconstructing a 3D Image of an object from cone 
beam projection data includes means, such as programmed computer, fordetenrtining values representing pla- 
10 nar integrals on a set of planes <^ in Radon space employing the procedures summarized above, and means, 
such as the pnDcessing apparatus disclosed in application Serial No. [RD-19564] for per- 

fomning an inverse Radon transfonn on the values representing the planar integrals on the set of planes to 
reconstruct an image of the object 

In the more particular embodiments where the planes comprise a set of coaxial planes containing a refer- 
15 ence axis intersecting the origin, the step of determining values representing planar integrals on the set of 
planes preferably includes the nested steps of, for each source position Si not on the reference axis, 

defining in Radon space a corresponding circle G/on the conresponding spherical shell in a plane 
containing the particular source position Sf and perpendicular to the planes intersections of the planes 
and the circles with the particular circle G/ defining on the circle G/ a plurality of points Pg corresponding to 
20 the circles D^, 

projecting the corresponding circle G/from the particular source position S/to a line Af,on the nor- 
malized detector plane, the points Pg projecting to conresponding points C^on the line M,, and 
for each projecting point Cg on the normalized detector plane, 

constructing lines 1/;^ on the normalized detector plane at a plurality of orientations passing 
25 through the projected point, the lines L/^^ being intersections on the nomialized detector plane of corresponding 
integration planes 0^,^ each containing a rotation axis along a line passing through the particular source position 
Sf, the particular point Pg, and the particular projected point Cg, 

rotating each of the lines on the normalized detector plane about the projected point 
Cg by a small angle 89 to define a line Lgj^ which Is the intersection of a plane Qg^' containing the particular 
30 rotation axis with the nomfiailzed detector plane, and detennining the rotation angle 8p between the planes 
and Qgk' by geometry from the angle 69, 

integrating along the lines LgK and L^' to detennine respective weighted line integrals 

and Jgic'f and 

dividing the difference between the weighted line Integrals Jgi^ and Jgi/ by the rotation angle 
35 6p to yield the radial derivative of the Radon datum at a point on the circle D,j where the plane Qgi, intersects 
the circle D^. . . 

In the more particular embodiments where the planes <l>j comprise a set of coaxial planes containing a refer- 
ence axis intersecting the origin, the step of determining values representing planar integrals on the set of 
planes <|)j includes the nested steps of, for each source position S| on the reference axis, 
40 for each plane <|>j intersecting the spherical shell corresponding to the particular source position 

Sf and defining a particular circle Dg, 

projecting the particular circle D/y from the particular source position S/ to a line Lg* on the 
normalized detector plane, 

constructing parallel lines Lgk on the normalized detector plane perpendicular to the line 
45 Lg*, the lines Lfp^ being Intersections on the normalized detector plane of corresponding integration planes 
each containing a rotation axis along a line passing through the particular source position S/and orthogonal to 
the plane of the particular circle Dg, 

translating each of the parallel lines Lip^ by a small distance to define a line L^' which is 
the intersection of a plane Q^^' containing the particular rotation axis with the normalized detector plane, and 
50 determining the rotation angle 8p between the planes and Q^^' by geometry from the distance between the 
lines Lgii and 

integrating along the lines Lg}f and Z.^' to detemnine respective weighted line integrals 

and Jgjt\ and 

dividing the difference between the weighted line Integrals and J^' by the rotation angle 
55 sp to yield the radial derivative of the Radon datum at a point on the drde Dg where the plane 0^^ intersects 
the circle Dg. 
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Brief Description of the Drawings - - — - 

While the novel features of the invention are set forth with particularity in the appended claims, the Inven- 
tion, both as to organization and content, will be better understood and appreciated from the following detailed 
5 description taicen in conjunction with the drawings, in which: 

FIG. 1 depicts a cone beann scanning geometry for 3D CT connected to reconstruction apparatus embody- 
ing the invention; 

FIGS. 2a, 2b. 2c, 2d, 2e and 2f are diagrams depicting the Radon transfonm approach to 3D CT imaging; 

FIG. 3 is a representation of the 3D Radon transform of an object at a given point; 
10 FIG. 4 depicts a set of coaxial planes <|>j in Radon space each containing a vertical or reference axis on 

which Radon data (planar Integrals) are to be detenmined; 

FIG. 5A depicts an object frame of reference and coordinate system; 

FIG. 5B depicts a normalized detector firame of reference and coordinate system; 

FIG. 5C depicts the manner in which the origins of the coordinate systems of FIGS. 5A and 5B coincide; 
IS FIG. 6 depicts an Integration plane, an integration frame of reference, and a conresponding coordinate sys- 
tem; 

FIG. 7 depicts cone beam data corresponding to an integration plane through the object; 
FIG. 8 similariy depicts cone beam data corresponding to a pair of closely spaced adjacent integration 
planes through the object; 
20 FIG. 9 shows geometry on an integration plane; 

FIG. 1 0 illustrates a procedure in accordance with the invention for computing the radial derivative of Radon 
data from x-ray cone beam data; 

FIG. 1 1 illustrates the Radon data on a plane that can be computed from one x-ray source position; 
FIG. 1 2 depicts a spherical shell or Radon shell representing all the Radon data that can be computed from 
25 one source position; 

FIG. 13 depicts an operation tenmed rotation about the a axis, as an example of Case 1; 
FIG. 14 depicts the Radon data that are generated by rotation about the a axis; 

FIG. 15 illustrates the generation of Radon data by performing the rotation operation on every point on a 
line projected onto the normalized detector plane; 
30 FIG. 1 6 depicts the rotation operation about the b axis, as an example of Case 2; - 

FIG. 17A depicts the Radon data that are generated by rotation about the 6 axis; ; - 

FIG. 178 depicts data points and lines on the nonmalized detector plane resulting from rotation about the 

b axis; 

FIG. 18A, line FIG. 4, depicts a set of coaxial vertical planes in Radon space; 
35 FIG. 18B depicts the Radon shell generated for a particular x-ray cone beam source position; : ~ 
FIG. 18C depicts the intersection of the FIG. 18A coaxial vertical planes with the Radon shell of FIG. 188; 
FIG. 19 represents in detail one of many circles where the vertical planes intersect the Radon shell as in . 
FIG. 1 8C, and illustrates the Case 1 procedure for generating Radon data on the vertical planes from detec- 
tor data; 

40 FIG. 20 represents the circles where the vertical planes intersect the Radon shell when the source position . 
S is on the axis of the vertical planes, and illustrates tfie Case 2 procedure for generating Radon data from 
the detector data; 

FIG. 21 depicts the intersection between two orthogonal planes. 
45 Detailed Description 

Referring initially to FIG. 1, a typical scanning and data acquisition configuration employing cone beam 
geometry connected to reconstruction apparatus embodying the invention. An object 20 is positioned within a 
field of view between a cone beam x-ray point source 22 and a 2D detector array 24, which provides cone beam 

so projection data. An axis of rotation 26 passes through the field of view and object 20. A midplane 28 may be 
defined which contains the x-ray point source 22 and is perpendicular to the axis of rotation 26. By convention, 
the axis of rotation 26 is referred to as the z-axis, and the intersection of the axis of rotation 26 and the midplane 
28 is taken as the origin of coordinates, x and y axes lie in the midplane 28 as indicated. For scanning the object 
20 at a plurality of source positions Si. the source 22 moves relative to the object 20 and the field of view along 

55 an appropriate scanning trajectory 30, while the detector 24 remains fixed with respect to the source 22. In FIG. 
1, the scanning trajectory 30 is for convenience of illustration shown as a circular scanning trajectory 30 lying 
the midplane 28, but other scanning trajectories may be employed and in fact are preferable, as is briefly dis- 
cussed hereinbelow. 
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The detector array 24 is connected to a data acquisition system (DAS) 32. During operation, x-ray photons 
that penetrate the object are detected by x-ray detector array 24 and registered by the data acquisition system 
(DAS) 32. The photon counts, after being nonmallzed by the air signals and converted to the negative of the 
logarithms, represent the line integrals through the object 20. Thus data are acquired at a number of source 

5 positions S, around the object 20 by scanning the source 22 and detector 24 along the scanning trajectory 30 
(or equivalently rotating the object 20 while the source 22 and detector 24 remain stationary). 

It should be noted, however, that the data set collected in such a single scan is incomplete, and artifacts 
may accordingly be introduced into the reconstruction, which may or may not be acceptable, depending on the 
particular application. Smith (1985, above) has shown that a cone beam data set is complete if there is a point 

10 from the x-ray source scanning trajectory on each plane passing through the object of interest (with the assump- 
tions that the detector is lodged in position relative to the source and large enough to span the object under 
inspection). A configuration suggested by MInerbo (1979. above) and by Heang K. Tuy, "An Inversion Formula 
for Cone-Beam Reconstruction". SIAM J. Math.. Vol. 43. No. 3, pp. 546-552 (June 1983), which Smith points 
out satisfies his condition for data completeness, is to employ two circular source scanning trajectories which 

15 are perpendicular to each other. Another scanning configuration which achieves data completeness is dis- 
closed in commonly-assigned U.S. Patent application Serial No. 07/572,651. filed August 27. 1990. by Eberhard 
et al. and entitled "SQUARE WAVE CONE BEAM SCANNING TRAJECTORY FOR DATA COMPLETENESS 
IN THREE-DIMENSIONAL COMPUTERIZED TOMOGRAPHr. Alternatively, instead of acquiring a complete 
cone beam x-ray data set, the invention of the above-Incorporated application Serial No. 

20 [RD-19695] may be employed, using optically-acquired object boundary infonnation to iteratively con-ect for 
missing data during the inverse Radon transfonm process. 

The data acquisition system (DAS) 32 is connected to a representative processor 34 which serves to recon- 
struct a 3D image of the object 20 by calculating planar integrals on a set of planes from the line integrals through 
the object 20 In accordance with the invention, and performing an inverse Radon transform on the planar integ- 

25 rate to reconstmct a three-dimensional image of the object 20. A suitable image display 36 is connected to the 
representative processor 34 as an output device. 

Refemng now to FIGS. 2A through 2f and FIG. 3, represented in general is the Radon transform approach 
to 3D imaging. 

Specifically, the object Itself is defined in terms of its x-ray attenuation coefficient /(x,y,z) (FIG. 2a). The 
30 measured cone beam projection data then corresponds to a line integral of this function over the rad ia! direction 
X(^) = l f (r, 9, Zo)dr (FIG. 2b). The line Integrals of the detector data (also known as detector Integrals) are 
given by J X(e) c/G = j; / (r. 0, Zo) drdQ (FIG. 2C). In the parallel beam case, these detector integrals are simply 
equal to the Radon transform of the object In the cone beam case, however, the Radon transform is given 
instead by // /(r, 8, Zo)rdrdQ (FIG. 2d). 
35 The additional factor of r in the Radon transfomi integral results from the Jacobian of the coordinate trans- 
fonmation from Cartesian to polar coordinates. As depicted in FIGS. 2e and 2f. an Inverse Radon transfonn 
procedure reconsfaucts a 3D CT image from the detector integrals. Since direct inverse Radon transformation 
requires planar integrals of the object as input, a preliminary step of converting line integrals (cone beam detec- 
tor integrals) to planar integrals (Radon data) is required, to which the present invention is directed. 
40 As depicted in FIG. 3. the 3D Radon transfonn of an object at a point x^, y^, is given by the area integral 
of the x-ray attenuation coefficient over the plane passing through Xo. /o, Zo that is perpendicular to the line 
from the origin to Xo, yo, Zo, and can be expressed as 



45 ^U.-^o.^cO^ jjnx,y,z)da 



(1) 



For a 2D radon transfonn, the situation is similar, except that the integral is over a line, not over a plane. 
The planar integral can also be expressed as 
^ ^ R{s,n) = !d^rS(s-r;h)f{!) (2) 

where n = {sinQcos^, sin 9 sin^, cos^ is a direction vector characterizing the nomial to the plane; s is the distance 
of the plane from the origin; and f{t) is the 3D object 

In words, R(s,n) represents the integrated density of the object on the plane whose normal is ft and which 
is at a distance s from the origin. The planar integral R{s,h) is also referred to as Radon data. 
55 The inverse Radon transformation by which a 3D object /(f) can be reconstructed from its planar integrals 
R can be expressed as 

W = ^ //JrfWcos B)ds^R(s,n)S(s -rti) (3) 
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-~ As described in detail In the above-incorporated application Serial No. - {Rd-19564). 

the inverse Radon transformation expressed in Equation (3) can be achieved through a two-step process. Step 
1 comprises 2D CT image reconstructions on a number of vertical planes (t>j in Radon space each containing 
the z axis. Step 2 comprises 2D CT image reconstructions on a number of horizontal planes. 
5 Thus, as represented in FIG. 4, what is required as input to the inverse Radon transfonmation are planar 

Integrals determined and organized on a plurality of planes <t)j containing a reference axis in Radon space, for 
example, on vertical planes 40, 42. 44 and 46 containing a vertical reference or z axis 48. 

Described in detail herein are methods and apparatus for converting the x-ray cone beam data to the planar 
integrals on the set of coaxial vertical planes or on any arbitrary set of planes, as input for inverse Radon 

10 transfomnation. 

With reference to FIGS. 5A. 5B and 5C, we first define the frames of reference and their associated coor- 
dinate systems used in the analysis herein. (It should be noted that the frames of reference and coordinate 
systems in the detailed discussion hereinbelow differ from the generalized coordinate system depicted in the 
representative scanning configuration of FIG. 1). In particular, FIG. 5A depicts an object frame of reference 

15 which Is fixed with respect to the object 20. Spatial location of any point In this frame of reference is expressed 
by the (x,y.z) triplet. The variables in Equations (2) and (3) are variables in the object frame of reference. Fig 
5B depicts a normalized detector frame of reference which is fixed with respect to a nonmallzed detector 50. 
Spatial locations in the normalized detector frame of reference are expressed by the (u,v,w) triplet. The origin 
of the (u,v,w) coordinate system is located at the center of the nonmalized detector 50, and the u and v axes 

20 lie on the normalized detector plane. The source position S\ is always on the w axis, but its distance from the 
center of the normalized detector 50 may vary from one source position S| to another source position Sj. 

As represented in FIG. 5C, we assume that the origin of the (x,y,z) coordinate system and the origin of the 
(u,v,w) coordinate system always coincide. In practice, this amounts to no more than scaling the actual detector 
24 readings to the plane that passes through the origin of the (x,y,z) coordinate system and is orthogonal to 

25 the line connecting the source position S/ and the origin. 

To facilitate the manipulation of planar integrals, we now introduce a third frame of reference. With refer- 
ence to FIG. 6, Q is any plane containing the source S. Let Plane Q intersect the nonmallzed detector plane at 
line L. Define a coordinate system (a,6,c) with source S as the origin ^such that d is a unit vector in plane Q 
orthogonal to line L, /) is a unit vector in plane Q parallel to L, and c = bxa. Whe shall refer to the coordinate 

30 system {a,b,c) as the coordinate system in the integration frame of reference. To further facilitate integration 
on plane Q, we note that each point on plane O is characterized by doublet {a,b) because the coondinate c is 
always zero. The doublet (a,i>) can be converted to polar coordinates (r.O) relative to the a,b axes by making 
the coordinate transformation: (a,6,c) (r,G,c), where 

r = Va2 + b^ 



35 
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-1^ 



9 = tan- 
a 



FIG. 7 illustrates a typical situation in cone beam scanning. Consider a plane Q, or slice, of the object as 
illustrated in FIG. 7. The cone beam source S projects plane Q on the normalized detector plane in the form 
of a straight line L In the Integration frame of reference, let the a axis intersect line L at point C (FIG. 6). By 
^ construction, SO is orthogonal to line L Let I SO I denote the distance between source S and the intersection 
point C. The datum X(t) on line t, where t represents the displacement from point C along L, is given by 



e-tan"' 



drde 



In other words, the datum X(i) represents the line integral of the object density along the line on plane Q 
making an angle 6 = tan-^(t/ 1 SO I) with the a axis. And noting that the variable / is proportional to tane, one 
would expect to obtain the integrated value of ffr.G.O) over the rand the 9 variables on plane Q, i.e., // /(r, 9, 
0)drdQ, by integrating X(/) (with some suitable weighting), over the / variable on line L. To this end let us express 
the quantity / = JJ /(r, 9, 0)drdQ in terms of the variable t. Now 

t = |sc|tan9 
dt = I SC I sec29d9 

Therefore the Integral / is given by 
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/ = jj f{r,e,0)drdd 



- f \SC\ , r 



(4) 



Thus the quantity / = // /(r. 9, 0)drrfe can be obtained by integrating the cone beam data X(0 on line L with 
weighting. In contrast the Radon data for this plane. R(s,n), where in the object frame of reference, s Is the 
distance of plane Qfitom the origin and n is its normal, is given by . 

^ = /i/(r.e.O)n/n/e (5) - 

Since / ^ R(s,n), the Radon data cannot be obtained by integrating cone beam data along straight lines on the 
normalized detector plane. (However, if the data on the detector were generated by parallel beams of x-rays, 
integrating data along straight lines on the normalized plane detector would yield the Radon data.) 

The only difference between the Radon data R{s,n) and the integral / is the absence of the factor r in the 
integral of /. Since the density value /(r.e.O) at each value of r is not known (otherwise there would be no need 
to do the cone beam scanning in the first place), the difference cannot be compensated for by weiqhtinq the 
data with r. ^ a 

One way to introduce the factor r into the integrand is to note that in rotating plane Q about any axis on 
the plane through the origin, each point on the plane is translated by an amount proportional to reiny where r 
IS the radial coordinate of the point, and y is its angular coordinate relative to the rotation axis. This observation 
prompts the following development with reference to FIGS. 8 and 9. 

Referring now to FIG. 8. let us consider another plane Q in the object which is very close to plane Q of 
FIG. 7. Plane Q- is obtained by rotating plane Q itself by a small angle 8p about a rotation axis a' on plane Q 
passing through the source position S. Selection of the actual rotation angle 8p is a compromise between accu- 
racy and signal-tonoise ratio in each particular system. A smaller 5p results in greater accuracy but at the 
expense of magnifying the noise in the data, and vice versa. Plane Q projects another line U on the normalized 
detector plane. Lines L and U interaect at point C. where the rotation axis sf intersects the normal detector 
plane. 

Let 59 be the angle between L and U. (It will be apparent that the angles 50 and 6p are related to each 
other by straightforward geometry, and that 5p can be detennined from 59. Examples of this calculation are 
provided hereinbelow for two particular cases.) Each point (r,9.0) on plane Q, in polar coordinates relative to 
the integration frame of reference, can be thought of as being obtained by translating by an amount 5rfrom the 
conresponding point in plane Q. Denoting by a the angle between the rotation axis a' and the a axis, the amount 
of translation Sr at point /^(r.O.O) Is given by 

Sr = S/3a'xr 

= 5/?(cos Ofl + sin aft) X (f cos Ba + r sin 0)6) 
= rsm{e~a)Spc 

The changes 8r in spatial coordinates r induce corresponding changes in the density values /(/), which in turn 
causes a change 8/ in the value of the integral /, which is gh^en by 
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81 = JJ V/(7. • Srdrdd 



Thus 



,5 ^ = £^^[/(^e.c),.„sin(0-a)]^d^(^e 

dc 

20 where 

T(c) = J/[/(r,e,c)sin(eKi)]rdrore 
The quantity 7(c) is almost the same as the Radon data in Equation (5) except for the extra factor sin(e - 
a). Now, however, the extra factor can be compensated for because it only involves the angular variable, whose 
value Is ava liable by measurement. This can be achieved by defining a new integral J wh ich includes a weighting 
25 factor in the integrand of / to cancel out the extra factor , 

J = rrlMm^ (6) 

•'•^ sin(0-o) ^ ' 
The geometry on plane Q Is shown in FIG. 9. From FIG. 9 we get 

sin e = . f 

cose =1^1 = 

- AC 
cosa = 4^ 



where AC denotes the displacement of C and C. Hence we have 
sin(8 - a) = sin 0COS a - cos 0 sin a 



(7) 



Now we can express the desired integral J in the variable t incorporating the weighting factor 1/sin(0 - a). 
Substituting Equations (4) and (7) into Equation (6) we have 

U fir,e,0)drd9 
sin(6^^a)^ 



55 i.^- ACJVPcf 



(r-AC)V|5C|' 



Going through the same mathematics as with integral / before, we obtain 
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dJ d rr 

ds 



(9) 



where, in the object frame of reference, s is the distance between plane Q and the origin, and 

'^(s.c) = fffirAOydrdQ 
10 is the planar integral of the function f on plane Q. 

Using Equation (9) we can in principle compute the radial derivative of the Radon data from the cone beam 
data and the Radon data themselves can be obtained by Integrating the result in the radial dimension. The pro- 
cedure Is illustrated In FIG. 10. To evaluate the radial derivative of the Radon data at a point P = sh where, in 
the object frame of reference, s = | OP I and n = OP I OP | , we do the following: 
15 1 . Determine the plane Q passing through the point P and orthogonal to line OP. 

2. Determine the line L where plane Q intersects the normalized detector plane. 

3. Locate the point ConL such that line SC is orthogonal to L 

4. Take any point C on line L, defining a rotation axis a' as a line from S to C 

Equivalently rotate plane Q about the rotation axis a' through a small angle 5p resulting in plane cy, and 
20 rotate line L about point C through a small angle 89 on the detector plane resulting in line L\ the plane Q' 
intersecting the nonnalized detector plane at the line L'. 

5. Compute the quantities J and J' on lines L and L\ respectively, using Equation (8). 

6. Compute the angle 6p from 69 by geometry. 

7. The radial derivative of the Radon data at point P is obtained from the quantities J, J, and 8p using the 
25 following equation: 

dR(s,n) ^ JW 
ds 6p 

Using the above procedure we can obtain the Radon data for all the planes through the object irradiated 
by the cone beam source. Incidentally, this is in agreement with the condition stated in Smith (1989, above) 

30 that. In order to have complete data In cone beam scanning every plane through the object should intersect a- 

source position. , ^ ^ , _ _ 

With reference now to FIGS. 11, 12 and 13. the range of the Radon data (i.e., planar Integrals) that can 
be generated at one source position in this way can be estimated quantitatively. Let the plane of FIG. 11 be 
any plane containing the cone beam source S and the origin O; call this plane W. Consider any plane Q ortho- 

35 gonal to plane IV and containing the source S. Let U be the line where plane Q intersects plane W; that is. 
plane Q passes through line U and Is orthogonal to plane the plane of FIG. 11. Let V be the line on plane 
ly passing through the origin and orthogonal to line L/, and let lines U and V intersect at point P. Let v be the 
unit vector along line V. As shown below in Appendix A with reference to FIG. 21 . the vector I OP | v Is orthogonal 
to plane Q, and therefore the planar integral over Q is the Radon datum P( I OP | , v) in the object frame of refer- 

40 ence, Le., a Radon datum at point P. Since the angle OPS Is a right angle, point P lies on the circle on plane 
ly with OS as diameter. By applying the same operation to all the planes orthogonal to plane IV and passing 
through source S, Radon data are generated on the entire circle, as Illustrated in FIG. 11. 

As represented in FIG. 12. by repeating the entire operation perfonned on plane W on all the other planes 
containing the line segment OS, Radon data are generated on all the circles containing OS as diameter. In other 

45 words. Radon data are generated on a spherical shell with OS as diameter. This spherical shell may also be 
refen-ed to as the Radon shell. 

Two particular cases will now be considered in detail, distinguished by the orientation of the rotation axis 
about which the FIG. 6 integration plane Q is rotated. As described hereinabove, each of the many rotation 
axes passes through the source position S. the line SO passing through the source position S and the origin 

50 O is orthogonal to the nomrtalized detector plane, and the nonnalized detector plane contains the origin. 

Case 1 applies when the rotation axis is not coincident with the b axis, and includes the specific case where 
the rotation axis is coincident with the a axis. Thus Case 1 may be described as extended rotation about the 
a axis. In the particular embodiments described herein where Radon data are being generated on a set of 
coaxial planes <1)| each containing a vertical or reference axis, as depicted in FIG. 4, Case 1 applies for all source 

55 positions S/ not on the vertical or reference axis. 

Case 2 applies when the rotation axis Is coincident with the b axis of FIG. 6. In this case, the rotation axis 
Is parallel to the normalized detector plane. In the particular embodiments described herein where Radon data 
are being generated on a set of coaxial planes ^ each containing a vertical or reference axis as depicted in 
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FIG. 4, Case 2 applies for all source position S/ which are on the vertical or reference axis. - 
A specific Case 1 case will not be considered, where the rotation axis is coincident with the a axis. In this 
case a = 0 in FIG. 10, and the two lines L and U Intersect at point C where the a axis Intersects the detector 
plane. If lines are drawn on the nonmalized detector plane at ail orientations passing through point C, (FIG. 1 3), 
5 they are the projections on the normalized detector plane of integration planes at different orientations contain- 
ing line SC as the common axis. Label these integration planes Q^^OzfOs , etc. Then according to Equation 

(9) the quantity dJ/dp, with the weighting function sine in J computed for each pair of adjacent lines closely 
spaced on the detector plane, yields the derivative of the planar integral on plane Q/ which projects onto the 
pair of adjacent lines. 

10 This situation is illustrated in FIG. 14, which depicts the view in direction SC, i.e., from the source toward 
the intersection point Con the detector. The plane of FIG. 14 contains the origin O, and line SC is orthogonal 
to FIG. 14. Point Pis the intersection point between the plane of FIG. 14 and line SC. Because line SC is ortho- 
gonal to the plane of FIG. 14, all the planes Q/ that contain line SC appear as lines forming the polar grid with 
point P as center; the lines are labeled A^Ai^A^,..., etc. in FIG. 14, con^esponding to planes Oi, 02,03 etc., 

15 respectively. From the origin O drop orthogonal lines to each of these lines, and let each pair of orthogonal 

lines intersect at locations 8^,82,8^ as illustrated. Again, as shown in Appendix A, the lines from the origin 

orthogonal to lines A^, 42,^3,.... are also orthogonal to planes Qi ,02,03 etc. Therefore the planar integrals 

on planes Qi, 02,03 etc. comprise the Radon data on points 81,82,8^ etc. And, because each 8/ Is at the 

intersection of two orthogonal lines, one from origin O and one from P, all points B/ fall on the circle with line 

20 segment OP as diameter. 

Since the points 8, all lie on the plane of FIG. 14, the plane of the circle on which they fall is orthogonal to 
line SC. Furthenmore, because line segment OP is orthogonal to line segment SP, point P lies on the surface 
of the sphere with OS as diameter. Thus P is the point where line segment SC intersects the Radon shell. 
The operations indicated in FIG. 13 may be summarized: 

25 (1) construct lines on the normal detector plane at all orientations passing through a point C, 

(2) compute the quantity J with the weighting function sinG on each of the lines, and 

(3) compute the derivative d J/dp. 

As a result, Radon data are generated on a circle (FIG. 14) on the plane containing the origin O and ortho- 
gonal to line SC, with line segment OP as a diameter of the circle where P is the point where line SC intersects 
30- the Radon shell. This entire operation is referred to herein as the rotation operation at point C on the detector 
plane. 

FIG. 15 represents the nonnalized detector plane for a particular source position. The rotation axis 
intersects the normalized detector plane at point C/, and lines 62, 64, 66, 68, 70 and 72 are representative 
intersections of various integration planes O with the normalized detector plane. To illustrate the rotation oper- 

35 ation in the context of FIG. 1 5, take any line M on the normalized detector plane and perform the rotation oper- 
ation on each point Cj on the line. For each point Cj a circle Dj of Radon data is generated, where the plane of 
the circle is orthogonal to line SC; and the diameter of the circle is line segment OPy, where Pj is the point where 
line SCy intersects the Radon shell. Since alt the points Cj line on line M, all points Pj lie on circle G where the 
plane containing source S and line M intersects the Radon shell. Also, because the plane of each circle Oy is 

40 orthogonal to the conresponding line SPy which lies on the plane of circle G, the plane of circle Oy is orthogonal 
to the plane of the G. 

To summarize, if rotation operations are performed on all the points on line M on the normalized detector. 
Radon data on a series of circles Dy are generated. The plane of circle Dy is orthogonal to the plane of circle G 
where the plane containing source S and line M intersects the Radon shell, and the diameter of Dj is line seg- 
45 ment OPy, where Py is a point on circle G. If points Cy are sufficiently finely sampled on line M, the set of Radon 
circles Dy generated is sufficiently dose to cover the entire Radon shell. In other words, the data on the entire 
Radon shell can be generated by performing rotation operations on a line on the detector plane. 

The weighting function sin9 used in computing the function J in this case is singular at G = 0. The singularity 
can be removed by the process of regularization, which has been successfully employed in the filtering portion 
50 of filtered backprojection in conventional x-ray CT where the kernel also contains singularity. 

Case 2 will now be considered where the rotation axis is coincident with the b axis in FIG. 6. In this case 
a = 7t/2 In FIG. 10, and the two lines L and U are parallel to each other since the rotation axis is parallel to the 
normalized detector plane. 

Refemng to FIG. 16, if all the locations along the line through DC, which is orthogonal to L (Appendix B), 
55 lines are drawn parallel to L and L', they are the projections on the detector plane of planes at different orien- 
tations containing the b axis as a comnK)n axis. Label these planes 0|, 02,03,..., etc. Then according to Equation 
(9) the quantity dJ/dp, with the weighting function cosG in J computed for each pair of adjacent lines closely 
spaced on the normalized detector plane, yields the derivative of the planar integral on the plane Oy which pro- 
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- jects onto the pair of adjacent lines. — _ 

This situation is illustrated in FIG. 17A. which depicts the view in the direction of the b axis, and in FIG 
17B which depicts data points and lines on the nonnalized detector plane for Case 2. The plane of the FIG 
17A is ttie c - a plane containing the source; in Appendix B It is shown that the origin 0 also lies on this plane' 
Since the b axis is orthogonal to the plane of FIG. 17A. all planes Qj that contain the b axis as a common axis 
appear as the lines fbmiing the polar grid with the source S as center; the lines are latieled as /\, A2A3 etc 

in FIG. 17A. oon^sponding to planes Q^.Q^.Q^ etc.. respectively. From the origin drop orthogonal lines to 

each of these lines, and let each pair of orthogonal lines intersect at locations 81.82,83 as illustrated Again 

as shovm in Appendix A, the lines from the origin orthogonal to the lines In the polar grid, A2M are also 

orthogonal to planes QuQ2.(h etc. Therefore, the planar integrals on Qu<h.Qz constitute the'ltedon data 

on points 81.82,6, eta And. because each Bj is at the intersection of two orthogonal lines, one from origin 

0 and one from source S. these points fall on the circle with line segments OS as diameter 

Since points Bj all lie on the plane of FIG. 1 7A. the plane of the cirde of Radon data is orthogonal to the b 
axis which IS nomial to the plane. Now, as shown In FIG. 178. the b axis Is parallel to the set of parallel lines 
on the detector plane including L and f. and the plane of the circle is orthogonal to the set of parallel lines. 

The operations indicated in FIG. 16 may be summarized: 

(1) construct lines parallel to direction 9 on the detector plane at all locations covering the enUre plane, 

(2) compute the quantity J with the weighting function cos9 on these lines, and 

(3) compute the derivative d;7dp. 

As a result. Radon data are generated on a circle (FIG. 17A) on the plane containing origin O and source S 
and orthogonal to the set of parallel lines on the detector plane, with line segments OS as a diameter of the 
circle. This entire operation is refenred to herein as the translation operation at angle 9. 

To generate Radon data on the entire Radon shell, translation operations are perfomied at all angles on 
the normalized detector plane. For each angle 9j a circle D, of Radon data is generated, where the plane of the 
circle is orthogonal to the lines at angle 6, on the detector plane and contains origin O and source S and the 
diameter of the circle is line segment OS. If the angles are sufficiently finely sampled, the set of Radon circles 
generated are sufficiently close to cover the entire Radon shell. 

Finally we can tackle the task we set out to resolve: generating Radon data on the set of coaxial vertical 
planes 9j of FIG. 4 from the cone beam data. 

In general, the procedure involves nested steps, designated hereinbeiow by subscripts I, j and k. Subscript 

1 conresponds to the various source positions S,. Thus, the subsequent steps are repeated for each of the source 
positions S, Subscript j corresponds to the set of planes 9j. which may be an arbitrary set of planes, in Radon 
space on which it is desired to generate planar integrals. For each of the source positions S, a corresponding 
Radon shell is defined on which Radon data can be detemiined, and inteisectlons of the planes e, with the 
Radon shell define a set of circles D^on the Radon shell. For each particular source position S,the subscript 

- 1 remains fixed, while the subscript j varies to designate each drcle of the set of cirdes D,. 

Then, for each of the cirdes Dj. each of the further subsequent steps is repeated. Specmcally. a rotation 
axis IS defined as a line through the particular source positfon S, interseding the particular circle D„, and per- 
pendicular to the plane of the cirde Dj(and perpendicular to the corresponding plane On each of tlie rotation 
axes a set of integration planes is defined, the integration planes Q,,belng coaxial with the particular rotation 
axis and intersection of the particular cirde Oj to define the location of a Radon datum point /?s,forthat particular 
integration plane For each particular source position S, and drde D, the subscripts i and j remain fixed 
while the subscript k varies to designate each Integration plane It will be appredated that the overall pro^ 
oedure of nested steps involves a multiplicity of Individual integration planes Qj», each coiresponding to an indi- 
vidual Radon datum point 

Continuing with the general procedure, each of the multiplldty of integration planes Qn, Intersects the nor- 
malized detector plane (the orientation of which con-esponds tothe particular source position S,) on respective 
lines V Then, for each of the lines on the normalized detector plane, the corresponding integration plane 
Qj» IS rotated about the rotation axis by a small rotation angle to define a plane Q,/ intersecting the normalized 
detector plane on a corresponding line V- Equivalently. the line V is rotated about a point on the line (Case 
1) or translated to a parallel line (Case 2). depending on the particular source position S, to define the line W 
and the con-esponding integration plane Q^'. 

Finally, to detennine the radial derivative of the Radon datum at the particular point weighted line integ- 
rals Jj» and V are determined by integrating respectively along the lines % and V. and the difference be- 
tween the weighted line integrals is divided by the rotation angle sp. 

ConsWering now the specific situation where it is desired to generate planar integrals on a set of coaxial 
vertical planes as input to the Radon inversion procedure of the above-incorporated application Serial No 

[RD-19564], FIG, 18A. like FIG. 4, illustrates the set of coaxial vertical planes conteining 
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— the z axis as the common axis, where the vertical axis is taken to be the z axis in the object frame of reference. - 
FIG. 18B shows a general cone beam scanning situation, with the source S, the origin O, and the Radon shell 
generated at the source position, since the source position S Is not on the z axis, FIG. 18B is an example of 
Case 1 as defined hereinabove. FIG. 1 8C Illustrates the geometry, showing Intersections of the FIG. 1 8A planes 
5 <|>j with the FIG. 1 8B Radon shell for a particular source position S/. It can be seen that these intersections define 
circles on the Radon shell of points on the individual planes 

. FIG. 19 illustrates in particular how to generate Radon data on the vertical planes (|>j in the Case 1 situation 
of FIG. 1 8C, where cone beam data Is acquired at each of a plural Ity of source positions S/, each source position 
S/ resulting in a corresponding spherical shell or Radon shell, such as the Radon shell of FIGS. 188. 18C and 
fO 19. Also, for each source position S/, a corresponding circle G/ is defined on the corresponding Radon shell in 
a plane containing the source position S/and perpendicular to the plane 4>j> i.e., a horizontal plane orthogonal 
to the vertical axis. 

It can be shown that, for each source position S/, each of the vertical planes <|)j Intersects the Radon shell 

corresponding to the particular source position S/ in a circle which may be designated circle Diy The circles D/; 
intersect the circle G/ at points corresponding to the circles D^. (Stated alternatively, the planes intersect 
the circle G/ at points P/y corresponding to the particular source position S/ and plane <|>j.) The circles D/y pass 
through the origin O, with lines from O to as diameters. In particular, 

G/ = The circle where the horizontal plane passing through source position S/ intersects the Radon shell 
^ D/y = The circle where the vertical plane (|>j containing the z axis intersects the Radon shell 
H/ = The horizontal equatorial plane through the center of the Radon shell 

Since origin O is on both vertical plane and the Radon shell, it is on intersection circle D/y. By symmetry, hori- 
zontal equatorial plane H/ through the center of the Radon shell must bisect circle D^, which is on a vertical 
plane <|)j. Therefore, the center of D^ lies on horizontal equatorial plane H/. Let OP^ be a diameter of circle D^, 

2g The distance between point and horizontal equatorial plane H/ Is the same as that between origin O and H/. 
It follows that point P^ lies on circle G/, because the distance between the horizontal plane containing G; and 
plane H/ is the same as that between origin O and plane H/. 

Because it is desired to generate planar integrals on each vertical plane <|)j, and because, for each source 
position S/, only those planar integrals on the corresponding Radon shell can be generated, for each source 

^ position S/ the planar integrals on each of the circles D^ for the particular source position S/, S/ projects G/ onto 
the normalized detector plane as a corresponding line M/. The intersection points P^ project to corresponding 
points C/y on the line M/. 

Then, the rotation is performed on each of the points on the line M/. The rotation operation on each point 
C/y results in the generation of Radon data on the entire circle corresponding to the particular point C^ Thus, 
each circle D^in FIG. 19 is represented by the FIG. 14 circle having diameter OP, with point P in FIG 14 rep- 
resenting each of the intersection points P/y in FIG. 19. The rotation axis for each of the points C/y is a line through 
the particular source position S/, the corresponding point P^ and the point C^. When the rotation operation is 
performed for all points C/y on a particular line M/ corresponding to a particular source position S/, Radon data 
on the entire Radon shell for the particular source position S/ are generated, organized as the circles D^ on the 
^ planes <|)j. When the operations are repeated for each source position S/, Radon data Is filled In on all the desired 
vertical planes. 

Summarizing the rotation operation for each projected point C/y in FIG. 19, lines /./y/f are constructed on the 
normalized detector plane at a plurality of orientations passing through the projected point C^. These lines L/y/^ 
are intersections on the normalized detector plane of corresponding integrations planes as are represented 
^ in FIGS. 13 and 14. Each of the integration planes contains a rotation axis along a line passing through the 
particular source position S/, the particular point P/y, and the particular projected point C/y. A multiplicity of figures 
like FIG. 14 can be drawn, the plane of each figure orthogonal to a particular rotation axis. 

Each of the lines L/y/^ is rotated on the normalized detector plane about the project point by a small angle 
SO to define a line Lp^ which is the intersection of a plane Q^/r' containing the particular rotation axis with the 
^ normalized detector plane. From the angle 50, the rotation angle 5p between the planes Q/y^r and Is deter- 
mined by geometry. Then respective weighted line integrals and J^/^' are determined by integrating along 
the lines L/y/f and L^^ in the manner described hereinabove. Finally, the difference between the weighted line 
integrals J^/^ and J^)l is divided by the rotation angle 6p to yield the radial derivative of the Radon datum at a 
point on the circle D^ where the plane 0/y/r intersects the circle D/y. 

The detennination of the rotation angle 5p between the two integration planes, given the angle 50 between 
the two detector lines can be accomplished by a number of geometric procedures. The following is a geomet- 
rically-derived fomiula for the Case 1 situation: 
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6p = cos(i)£os 1+tan2e — 
M+cos2f|yiitan2e 

where 

= azimuthal angle of plane <|>j with respect to SO 
Tl = polar angle of SO 

e = angle between the line L^p, and the reference line where the plane Intersects the normal ized detector 
plane 

The above Case 1 procedure cannot be applied when the source position S, Is on the z axis either directly 
,0 above or directly below the origin O, which are Case 2 situations as defined hereinabove. 

Referring in particular to FIG. 20, illustrated is a situation where the source position S, is directly below the 
origin O. This situation is also represented in FIGS. 1 6, 1 7A and 1 7B, where each axis of rotation passes through 
the source position S, parallel to the plane of the normalized detector. In FIG. 20, each of the vertical planes 
<|>j with the z axis as the common or reference axis intersects the Radon shell in a great circle Df, with OS, as 
,5 the diameter. Each circle has a corresponding axis of rotation for the Integration planes which Is perpen- 
dicular to the circle and Its corresponding plane 

For Case 2. each particular circle Dg is projected from the source position Si to a line Lg* on the normalized 
detector plane. Because circle contains O and S^ line is also the projection of the plane ^ onto the nor- 
malized detector plane. 

2Q To generate Radon data on the circle D^, parallel lines £.^are constructed on the normalized detector plane 
perpendicular to the line L,/. The tines L^are represented by the parallel lines on the nomiallzed detector plane 
in FIGS. 16 and 17B, and are the intersections on the detector plane of corresponding integration planes Q^^^. 
each containing a rotation axis along a line passing through the particular source position S,and orthogonal to 
the plane of the particular circle D^. 

25 Each of the parallel lines Lip, is translated by a small distance to define a line V which is the intersection 
of a plane Qgi,' containing the particular rotation axis with the nonnalized detector plane. The rotation angle 8p 
between the two planes Qg^ and Q,/ is determined by geometry from the distance between the lines Lg„ and 
Then respective weighted line integrals Jfj„ and J^^' are detemiined by Integrating along the lines Lgf, and 
V 'n the manner described hereinabove. Finally, the difference between the weighted line integrals J^i, and 

30 Jfjk' is divided by the rotation angle 8p to yield the radial derivative of the Radon datum at a point on the circle 
Off where the plane 0^,^ intersects the circle Dg. 

The detennination of the rotation angle 6p between the two integration planes can be accomplished by the 
following geometrlcally-drlved formula for the Case 2 situation: 

35 ^ TsoI^T7 

where 

I SO I = distance between S and O 
/ = distance of the line from O 
5y = translation distance 

^ Considering the number of computations required in embodiments of the invention, at each view angle the 
generation of the datum at one point in the Radon space requires computations on a line of data on the detector 
plane, which contains «N data points. Therefore, to generate the data on a circle on the Radon shell requires 
« A/ X A/ = /V2 computations. The number of computations required to generate data on N circles covering the 
Radon shell is therefore equal to » A/ x W2= A/3. Finally, the total amount of computation at all N view anqles is 
^ givenbyA/3xA/=A/*. 

A faster way to anive at the same estimate is the following: To generate the datum at one point In the Radon 
space requires N computations. Because there are « hfi points in the Radon space, the total amount of com- 
putations Is equal to » A^ x A/ = A^. 

While specific embodiments of the invention have been Illustrated and described herein, it is realized that 
modifications and changes will occur to those skilled in the art. It is therefore to be understood that the appended 
claims are intended to cover all such modifications and changes as fall within the true spirit and scope of the 
invention. 
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APPENDIX ft 

With reference to FIG. 21, it will be shown that 
vector pP\v of FIG. 11 is orthogonal to plane Q which 
intersects plane W at line U in FiG.* 11. "Define two planes 
and Q2 such that : 

« unit normal to Qi 
= unite normal to 



Let Qi and Q2 intersect at line L. Without loss of generality 
20 let the origin be on line L. 

Since "h-^, and h,lQ,, therefore n^eQ^. 
In a similar manner one can show that nx^Qj. 
Let be a unit vector along line L. 



rt3 € L C Q2 => n^JJij 

Hence the set {n^.n^.n^} forms an orthonormal basis of the 

space. Also, c?j is spanned by fu^ and /t, , and ' i s spanned by 
rt, and hy. 

Let £ be any point on Qi. Then P^X^^ + X^n^ for some* 
scalars X^ and A3. Let £1 = X.'n,-^ X^% be the point on Q2 clos- 
est to £ on Q2. Now 

For a fixed R, the minimum of \PPf occurs at X^' = 0 and 
A'jSsAj, i.e., 
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is a point on line i. Now 

Thus it is obvious that E£l is orthogonal to line L and 
orthogonal to plane Q^. 

Analysis of the geometry of FIG. 20 in the Case 2 

Situation: . 

SC is along a 
L is along b 

Thus SC is orthogonal to L, 
OS is along w 

L is on the {u,v) plane 

Thus OS is orthogonal to L 

Since 5C is orthogonal to L and OS is orthogonal to we 
. . .conclude that the plane containing lines SC and OS is 
orthogonal to line L. Since . is along the ^ axis, this 
Plane is the (c,a) plane through S. This plane contains the 
40 points, 0,S, and C. 
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1. A method for constnjcting a 3D image of an object from cone beam projection data, the cone beam pro- 
jection data being in the forni of line integrals through the object organized for each of a plurality of x-ray 
source positions S, as a 2D data set on a normalized detector plane containing an origin and perpendicular 
so to a line from each particular source position S, to the origin, said method comprising: 

determining values representing planar integrals on a set of planes ^ in Radon spaced by, for each 
of the source positions S/, 

defining in Radon space a con-esponding spherical shell on which Radon data can be determined, 
Intersections of the planes cfj with the spherical shell corresponding to the particular source position S/defr 
55 nlng a set of circles Dg on the spherical shell, and 

for each of the circles D^ 

defining a rotation axis as a line through the particular source position S/, intersecting 
the particular circle Dy. and orthogonal to the plane of the particular circle D*, 
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defining a set of coaxial integration planes 0^;^ each of the integration planes Q/;/fCon-- 
taining the particular rotation axis and intersecting the particular circle Dyto defrne the location of a Radon 
datum point R^^, and the Integration planes Op Intersecting the normalized detector plane on respective 
lines Lp, and 

for each of the lines Lijk on the nomnallzed detector plane, 

rotating the conresponding Integration plane 0^;^ by a small rotation angle 5p 
to define a plane Op', intersecting the nomialized detector plane on a corresponding line L^', 

integrating along the lines Lp and L^' to determine respective weighted line 

integrals Jp and J^', and 

dividing the difference between the weighted line integrals and Jp' by the 
rotation angle 5p to yield the radial derivative of the Radon datum at the particular point R^; and 

perfonning an Inverse Radon transfonm on the values representing planar Integrals on the set of 
planes to reconstruct an image of the object 

A method In accordance with Claim 1, wherein 

the planes 4>j comprise a set of coaxial planes containing a reference axis intersecting the origin; 

and 

the step of detemnining values representing planar integrals on the set of planes 4>j comprises, for 
each source position S/ not on the reference axis, 

defining In Radon space a corresponding circle on the corresponding spherical shell In a 
plane containing the particular source position S/ and perpendicular to the planes (t>j, intersections of the 
planes <|)j and the circles with the particular circle G/ defining on the circle G/ a plurality of points cor- 
responding to the circles D,j, 

projecting the corresponding circle G/ from the particular source position S/ to a line M/ on 
the normalized detector plane, the points P/y projecting to corresponding points on the line M/, and 
for each projected point C^on the normalized detector plane, 

constructing lines Lp on the normalized detector plane at a plurality of orientations 
passing through the projected point, the lines Lp being intersections on the normalized detector plane of 
corresponding integration planes Op each containing a rotation axis along a line passing through the par- 
ticular source position Sf, the particular point P/y, and the particular projected point Cy, 

rotating each of the lines Lp on the normalized detector plane about the projected 
point C/y by a small angle $0 to define a line Lp' which is the Intersection of a plane Op' containing the 
particular rotation axis with the normalized detector plane, and detemnining the rotation angle between 
the planes Op and Op by geometry from the angle 56, 

Integrating along the lines Lp and Lp' to detennine respective weighted fine integrals 

Jp and Jp', and 

dividing the difference between the weighted line Integrals Jp and Jp' by the rotation 
angle 5p to yield the radial derivative of the Radon datum at a point on the circle Dy where the plane Op 
intersects the circle D^. 

A method in accordance with Claim 1 , wherein 

the planes <|>j comprise a set of coaxial planes containing a reference axis intersecting the origin; 

and 

the step of detenmining values representing planar integrals on the set of planes <|»j comprises, for 
each source position S/ on the reference axis, 

for each plane <{)j intersecting the spherical shell conresponding to the particular source position S/ 
and defining a particular circle D^, 

projecting the particular circle Du from the particular source position S/ to a line L^* on the 
nonmalized detector plane, 

constructing parallel lines Lp on the nonnalized detector plane perpendicular to the line Ly*, 
the lines Z.p being intersections on the nonnalized detector plane of conresponding integration planes Op 
each containing a rotation axis along a line passing through the particular source position S/and orthogonal 
to the plane of the particular circle D^, 

translating each of the parallel lines Lp by a small distance to define a line Lp which is the 
intersection of a plane Op' containing the particular rotation axis with the nonmalized detector plane, and 
determining the rotation angle 5p between the planes Op and Op' by geometry from the distance between 
the lines Lp and Lp', 

Integrating along the lines Lp and Lp' to detennine respective weighted line integrals Jp and 
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Jjpt'i and - 



dividing the difference between the weighted line integrals and V by the rotation angle 
5p to yield the radial derivative of the Radon datum at a point on the circle where the plane Qu, intersects 
the circle *' 

4. A method in accordance with Qaim 2. wherein the step of detennining values representing planar Integrals 
on the set of planes comprises, for each source position S, on the reference axis: 

for each plane ^ intersecting the spherical shell corresponding to the particular source position S, 
and defining a particular circle Dg, 

projecting the particular circle from the particular source position S, to a line L,* on the 
normalized detector plane, 

constructing parallel lines Lji^on the nonmalized detector plane perpendicular to the line Z.,,*, 
the lines being intersections on the nomfiallzed detector plane of conrespondlng integration planes Q^^ 
each containing a rotation axis along a line passing through the particular source position S,and orthogonal 
to the plane of the particular circle Dg, 

translating each of the parallel lines by a small distance to define a line L^' which is the 
intersection of a plane Q^/ containing the particular rotation axis with the normalized detector plane and 
determining the rotation angle 5p between the planes Q^^^and Q^' by geometry from the distance between 
the lines Lyjt and Ly;t', 

integrating along the lines /.^and L^' to determine respective weighted line integrals Jn,, and 

J^', and 

dividing the difference between the weighted line integrals J^p, and Jg^' by the rotation angle 
6p to yield the radial derivative of the Radon datum at a point on the circle where the plane a* intersects 
the circle Dg, * 

5. Apparatus for reconstructing a 3D image of an object from cone beam projection data, the cone beam pro- 
jection data being in the form of a line integrals through the object organized for each of a plurality of x-ray 
source positions Sf as a 2D data set on a nonnalized detector plane containing an origin and perpendicular 
to a Irne from each particular source position S/ to the origin, said apparatus comprising: 

means for detemilning values representing planar integrals on a set of planes <(>. in Radon space 
by, for each of the source positions S/, 

defining in Radon space a con-esponding spherical shell on which Radon data can be determined, 
intersections of the planes with the spherical shell corresponding to the particular source position S,defh 
ning a set of circles Dg on the spherical shell, and 
for each of the circles Dg, 

defining a rotation axis as a line through the particular source position Sf, Intersecting 
the particular circle Dg, and orthogonal to the plane of the particular circle Dg, 

defining a set of coaxial integration planes Q/yj^eachof the integration planes Qjk con- 
taining the particular rotation axis and intersecting the particular circle Dg to define the location of a Radon 
datum point and the Integration planes Qfj^ intersecting the nonmalized detector plane on respective 
lines Lgjt, and 

for each of the lines L^/^on the normalized detector plane. 

rotating the conrespondlng integration plane by a small rotation angle 8p to define 
a plane Q^^', intersecting the normalized detector plane on a conrespondlng line V, 

integrating along the lines Lg^ and L^^' to detenmine respective weighted line integrals 

Jflatand J^'. and 

dividing the difference between the weighted line integrals Jgj, and V by the rotation 
angle 8p to yield the radial derivative of the Radon datum at the particular point Rgj^ and 

means for performing an inverse Radon transform on the values representing planar integrals on 
the set of planes to reconstruct an Inriage of the object 

6. Apparatus in accordance with Claim 5. wherein 

and '^'^"^ * ^"^^"^^ ® °^ P*®"®® containing a reference axis intereecting the origin; 

said means for detemilning values representing planar integrals on the set of planes is operable, 
for each source position S/ not on the reference axis, 

to define in Radon space a correspnnding circle G,on the con-esponding spherical shell in 
a plane containing the particular source position S, and perpendicular to the planes Intersections of the 
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planes <|)j and the circles Dgw'ith the particular circle G/ defining on the circle G/a plurality of points cor- 
responding to the circles D^, 

to project the corresponding circle G/ from the particular source position S/ to a line M/ on 
the normalized detector plane, the points projecting to corresponding points C^on the line and 
for each projected point C^on the normalized detector plane, 

to construct lines Lfji^on the nomnallzed detector plane at a plurality of orientations pas- 
sing through the projected point, the lines LjfK being intersections on the normalized detector plane of cor- 
responding integration planes each containing a rotation axis along a line passing through the 
particular source position Si, the particular point P^, and the particular projected point C^, 

to rotate each of the lines L^^ on the nomnafized detector plane about the projected 
point Cu by a small angle 69 to define a line Lfjn' which is the intersection of a plane Q^' containing the 
particular rotation axis with the nonmalized detector plane, and determining the rotation angle 5p between 
the planes and Qgi! by geometry from the angle 50, 

to integrate along the lines L^and to determine respecth/e weighted line integrals 

Jukand Jp', and 

to divide the difference between the weighted line Integrals Jg^ and Jyff' by the rotation 
angle 5p to yield the radial derivative of the Radon datum at a point on the circle Dy where the plane 
intersects the circle D/y. 

A method in accordance with Claim 5, wherein 

the planes i|>j comprise a set of coaxial planes containing a reference axis intersecting the origin; 

and 

said means for determining values representing planar integrals on the set of planes (|)j is operable, 
for each source position S; not on the reference axis, 

for each plane <t)j intersecting the spherical shell conresponding to the particular source position S/ 
and defining a particular circle D/y, 

to project the particular circle from the particular source position 5/ to a line Ly* on the 
normalized detector plane, 

to construct parallel lines L,y;f on the nomrialized detector plane perpendicular to the line L,f, 
the lines Ly^ being intersections on the normalized detector plane of corresponding integration planes Qij^ 
each containing a rotation axis along a line passing through the particular source position S/and orthogonal 
to the plane of the particular circle D^, 

to translate each of the parallel lines Lgk by a small distance to define a line L^' which Is the 
intersection of a plane 0/;;^' containing the particular rotation axis with the normalized detector plane, and 
determining the rotation angle 5p between the planes 0/y/cand O// by geometry from the distance between 
the lines Z.^^ and L/y/, 

to integrate along the lines Z.^, and Lgi/ to determine respective weighted line integrals 

and Jul/, and 

to divide the difference between the weighted line Integrals J/yif and Jfjk by the rotation angle 
$P to yield the radial derivative of the Radon datum at a point on the circle where the plane intersects 
the circle Dg. 

Apparatus in accordance with Claim 6, wherein said means for detenmining values representing planar 
Integrals on the set of planes <t)j is operable, for each source position S/ on the reference axis: 

for each plane intersecting the spherical shell conresponding to the particular source position S/ 
and defining a particular circle Dy, 

to project the particular circle Dy from the particular source position S/ to a line Ly* on the 
normalized detector plane, 

to construct parallel lines L^j^on the normalized detector plane perpendicular to the line Ly*, 
the lines Ly^ being intersections on the normalized detector plane of corresponding integration planes 
each containing a rotation axis along a line passing through the particular source position S/and orthogonal 
to the plane of the particular circle Dy, 

to translate each of the parallel lines Lyi^ by a small distance to define a line Lyk' which Is the 
intersection of a plane Qy}/ containing the particular rotation axis with the normalized detector plane, and 
determining the rotation angle 5p between the planes QyK and Qyn' by geometry from the distance between 
the lines Lyi^ and Ly^, 

to integrate along the lines LyK and L/ji/ to determine respective weighted line integrals Jy^ 

and Jyi(', and 
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FIG. 3 
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FIG. 12 
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FIG. 13 




FIG. 14 
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